library(MendelianRandomization)
data_dir = "../../dat"
seq_predictions_all = read.csv(file.path(data_dir, "means_and_uncertainties_hopefully_correct.csv"))
seq_predictions = subset(seq_predictions_all, seq_num == 25)
str(seq_predictions)
'data.frame': 3001 obs. of 6 variables:
$ seq_num : int 25 25 25 25 25 25 25 25 25 25 ...
$ mut : Factor w/ 3000 levels "0A","0C","0G",..: 12 1 2 3 342 343 344 674 675 676 ...
$ X_pred_mean: num 5.66e-06 2.34e-05 4.55e-05 -3.47e-05 1.45e-04 ...
$ X_pred_var : num 3.17e-05 3.95e-05 6.99e-05 3.93e-05 9.92e-05 ...
$ Y_pred_mean: num 4.68e-06 1.52e-05 8.47e-05 -9.28e-06 8.67e-05 ...
$ Y_pred_var : num 2.52e-05 3.08e-05 8.28e-05 2.55e-05 9.15e-05 ...
bx <- unlist(seq_predictions["X_pred_mean"])
bxse <- unlist(seq_predictions["X_pred_var"])
by <- unlist(seq_predictions["Y_pred_mean"])
byse <- unlist(seq_predictions["Y_pred_var"])
MRInputObjectEx <- mr_input(bx = bx,
bxse = bxse,
by = by,
byse = byse)
cis = c()
for (seq in (1:50)) {
seq_predictions = subset(seq_predictions_all, seq_num == seq)
bx <- unlist(seq_predictions["X_pred_mean"])
bxse <- unlist(seq_predictions["X_pred_var"])
by <- unlist(seq_predictions["Y_pred_mean"])
byse <- unlist(seq_predictions["Y_pred_var"])
MRInputObject <- mr_input(bx = bx,
bxse = bxse,
by = by,
byse = byse)
EggerObject <- mr_egger(MRInputObject, robust = FALSE, alpha = .05)
# cis = append(cis, c(EggerObject$CILower.Est, EggerObject$CIUpper.Est))
print(seq)
print(c(EggerObject$CILower.Est, EggerObject$CIUpper.Est))
print(c(EggerObject$I.sq))
}
[1] 1
[1] 0.05687267 0.06067051
[1] 0
[1] 2
[1] 0.3306068 0.3624528
[1] 0.2105841
[1] 3
[1] 0.03796090 0.04176579
[1] 0.3214789
[1] 4
[1] 0.1549294 0.1735800
[1] 0.2925369
[1] 5
[1] 0.1533373 0.1661913
[1] 0.07879359
[1] 6
[1] 0.1685509 0.1815121
[1] 0
[1] 7
[1] 0.2675980 0.2861198
[1] 0
[1] 8
[1] 0.2909661 0.3147951
[1] 0.1831951
[1] 9
[1] 0.1449050 0.1556574
[1] 0.1622787
[1] 10
[1] 0.1257318 0.1356557
[1] 0
[1] 11
[1] 0.1473885 0.1603717
[1] 0
[1] 12
[1] 0.04852230 0.05193937
[1] 0
[1] 13
[1] 0.2326184 0.2606935
[1] 0
[1] 14
[1] 0.4096944 0.4500918
[1] 0
[1] 15
[1] 0.1006563 0.1091880
[1] 0.02123769
[1] 16
[1] 0.03630137 0.03925981
[1] 0.2420155
[1] 17
[1] 0.08023490 0.08748782
[1] 0.1503298
[1] 18
[1] 0.2139596 0.2311665
[1] 0
[1] 19
[1] 0.2752876 0.2981095
[1] 0
[1] 20
[1] 0.2259239 0.2427902
[1] 0
[1] 21
[1] 0.01434114 0.01553469
[1] 0
[1] 22
[1] 0.1928645 0.2083450
[1] 0
[1] 23
[1] 0.05310580 0.05783676
[1] 0
[1] 24
[1] 0.3177343 0.3388181
[1] 0
[1] 25
[1] 0.2673437 0.2889935
[1] 0.2421144
[1] 26
[1] 0.02250383 0.02425483
[1] 0
[1] 27
[1] 0.05180822 0.05611104
[1] 0
[1] 28
[1] 0.05388898 0.05843682
[1] 0
[1] 29
[1] 0.04815811 0.05224462
[1] 0
[1] 30
[1] 0.08851500 0.09372713
[1] 0.174372
[1] 31
[1] 0.3338592 0.3653642
[1] 0
[1] 32
[1] 0.03867074 0.04214653
[1] 0
[1] 33
[1] 0.4801847 0.5207160
[1] 0
[1] 34
[1] 0.03362065 0.03616059
[1] 0
[1] 35
[1] 0.05669501 0.06187179
[1] 0
[1] 36
[1] 0.1737785 0.1876079
[1] 0.1640396
[1] 37
[1] 0.08575466 0.09214800
[1] 0
[1] 38
[1] 0.06662938 0.07251736
[1] 0
[1] 39
[1] 0.1864075 0.2002264
[1] 0
[1] 40
[1] 0.2280670 0.2459672
[1] 0
[1] 41
[1] 0.02417243 0.02568200
[1] 0
[1] 42
[1] 0.08946369 0.09645143
[1] 0
[1] 43
[1] 0.1533693 0.1655198
[1] 0
[1] 44
[1] 0.04358466 0.04715561
[1] 0
[1] 45
[1] 0.1027557 0.1105777
[1] 0
[1] 46
[1] 0.1188304 0.1293053
[1] 0
[1] 47
[1] 0.02712174 0.02866630
[1] 0.05283552
[1] 48
[1] 0.1388547 0.1522203
[1] 0
[1] 49
[1] 0.4762848 0.5147806
[1] 0
[1] 50
[1] 0.09995734 0.10829652
[1] 0.09234253
# dim(cis) <- c(2, 50)
# print(cis)
mr_plot(MRInputObjectEx)
install.packages("devtools")
Installing package into ‘/home/stephenmalina/R/x86_64-pc-linux-gnu-library/3.6’
(as ‘lib’ is unspecified)
also installing the dependency ‘remotes’
trying URL 'https://cloud.r-project.org/src/contrib/remotes_2.1.1.tar.gz'
Content type 'application/x-gzip' length 134042 bytes (130 KB)
==================================================
downloaded 130 KB
trying URL 'https://cloud.r-project.org/src/contrib/devtools_2.2.2.tar.gz'
Content type 'application/x-gzip' length 375464 bytes (366 KB)
==================================================
downloaded 366 KB
* installing *source* package ‘remotes’ ...
** package ‘remotes’ successfully unpacked and MD5 sums checked
** using staged installation
** R
** inst
** byte-compile and prepare package for lazy loading
** help
*** installing help indices
** building package indices
** installing vignettes
** testing if installed package can be loaded from temporary location
** testing if installed package can be loaded from final location
** testing if installed package keeps a record of temporary installation path
* DONE (remotes)
* installing *source* package ‘devtools’ ...
** package ‘devtools’ successfully unpacked and MD5 sums checked
** using staged installation
** R
** inst
** byte-compile and prepare package for lazy loading
** help
*** installing help indices
*** copying figures
** building package indices
** installing vignettes
** testing if installed package can be loaded from temporary location
** testing if installed package can be loaded from final location
** testing if installed package keeps a record of temporary installation path
* DONE (devtools)
The downloaded source packages are in
‘/tmp/RtmpfKGeTB/downloaded_packages’
library(devtools)
Loading required package: usethis
install_github("WSpiller/RadialMR")
Skipping install of 'RadialMR' from a github remote, the SHA1 (18a7a6e9) has not changed since last install.
Use `force = TRUE` to force installation
plot(by/bx, bx/byse, xlim=c(min((by-2*byse)/bx), max((bx+2*byse)/bx)), ylim=c(0, max(bx/byse)))
for (j in 1:length(bx)) {
lines(c((by[j]-1.96*byse[j])/bx[j], (by[j]+1.96*byse[j])/bx[j]),
c(bx[j]/byse[j], bx[j]/byse[j]))
}
abline(h=0, lwd=1); abline(v=0, lwd=1)

library(RadialMR)
r_input <- format_radial(bx, by, bxse, byse)
funnel_radial(egger_radial(r_input, alpha = .05))
funnel_radial(ivw_radial(r_input, alpha = .05))
Weights not specified: Adopting modified second-order weights
Radial IVW
Residual standard error: 0.948 on 3000 degrees of freedom
F-statistic: 2820.08 on 1 and 3000 DF, p-value: <2e-16
Q-Statistic for heterogeneity: 2698.411 on 3000 DF , p-value: 0.9999713
Outliers detected
Number of iterations = 6
Error: Tibble columns must have consistent lengths, only values of length one are recycled:
* Length 4: Columns `x`, `xend`
* Length 3001: Column `PANEL`

plot_radial(egger_radial(r_input, alpha = .05))
LS0tCnRpdGxlOiAiTWVuZGVsaWFuIFJhbmRvbWl6YXRpb24gYW5hbHlzaXMgb2YgRGVlcFNFQSBnZW5lcmF0ZWQgZGF0YSIKb3V0cHV0OiBodG1sX25vdGVib29rCi0tLQpgYGB7cn0KbGlicmFyeShNZW5kZWxpYW5SYW5kb21pemF0aW9uKQpgYGAKCmBgYHtyfQpkYXRhX2RpciA9ICIuLi8uLi9kYXQiCnNlcV9wcmVkaWN0aW9uc19hbGwgPSByZWFkLmNzdihmaWxlLnBhdGgoZGF0YV9kaXIsICJtZWFuc19hbmRfdW5jZXJ0YWludGllc19ob3BlZnVsbHlfY29ycmVjdC5jc3YiKSkKc2VxX3ByZWRpY3Rpb25zID0gc3Vic2V0KHNlcV9wcmVkaWN0aW9uc19hbGwsIHNlcV9udW0gPT0gMjUpCmBgYAoKYGBge3J9CnN0cihzZXFfcHJlZGljdGlvbnMpCmBgYApgYGB7cn0KYnggPC0gdW5saXN0KHNlcV9wcmVkaWN0aW9uc1siWF9wcmVkX21lYW4iXSkKYnhzZSA8LSB1bmxpc3Qoc2VxX3ByZWRpY3Rpb25zWyJYX3ByZWRfdmFyIl0pCmJ5IDwtIHVubGlzdChzZXFfcHJlZGljdGlvbnNbIllfcHJlZF9tZWFuIl0pCmJ5c2UgPC0gdW5saXN0KHNlcV9wcmVkaWN0aW9uc1siWV9wcmVkX3ZhciJdKQpgYGAKCmBgYHtyfQpNUklucHV0T2JqZWN0RXggPC0gbXJfaW5wdXQoYnggPSBieCwgCiAgICAgICAgICAgICAgICAgICAgICAgICAgYnhzZSA9IGJ4c2UsIAogICAgICAgICAgICAgICAgICAgICAgICAgIGJ5ID0gYnksIAogICAgICAgICAgICAgICAgICAgICAgICAgIGJ5c2UgPSBieXNlKQpgYGAKCgpgYGB7cn0KY2lzID0gYygpCmZvciAoc2VxIGluICgxOjUwKSkgewogIHNlcV9wcmVkaWN0aW9ucyA9IHN1YnNldChzZXFfcHJlZGljdGlvbnNfYWxsLCBzZXFfbnVtID09IHNlcSkKICAKICBieHQgPC0gdW5saXN0KHNlcV9wcmVkaWN0aW9uc1siWF9wcmVkX21lYW4iXSkKICBieHNldCA8LSB1bmxpc3Qoc2VxX3ByZWRpY3Rpb25zWyJYX3ByZWRfdmFyIl0pCiAgYnl0IDwtIHVubGlzdChzZXFfcHJlZGljdGlvbnNbIllfcHJlZF9tZWFuIl0pCiAgYnlzZXQgPC0gdW5saXN0KHNlcV9wcmVkaWN0aW9uc1siWV9wcmVkX3ZhciJdKQogIE1SSW5wdXRPYmplY3QgPC0gbXJfaW5wdXQoYnggPSBieHQsIAogICAgICAgICAgICAgICAgICAgICAgICAgIGJ4c2UgPSBieHNldCwgCiAgICAgICAgICAgICAgICAgICAgICAgICAgYnkgPSBieXQsIAogICAgICAgICAgICAgICAgICAgICAgICAgIGJ5c2UgPSBieXNldCkKICBFZ2dlck9iamVjdCA8LSBtcl9lZ2dlcihNUklucHV0T2JqZWN0LCByb2J1c3QgPSBGQUxTRSwgYWxwaGEgPSAuMDUpCiAgIyBjaXMgPSBhcHBlbmQoY2lzLCBjKEVnZ2VyT2JqZWN0JENJTG93ZXIuRXN0LCBFZ2dlck9iamVjdCRDSVVwcGVyLkVzdCkpCiAgcHJpbnQoc2VxKQogIHByaW50KGMoRWdnZXJPYmplY3QkQ0lMb3dlci5Fc3QsIEVnZ2VyT2JqZWN0JENJVXBwZXIuRXN0KSkKICBwcmludChjKEVnZ2VyT2JqZWN0JEkuc3EpKQp9CiMgZGltKGNpcykgPC0gYygyLCA1MCkKIyBwcmludChjaXMpCmBgYAoKCmBgYHtyfQptcl9wbG90KE1SSW5wdXRPYmplY3RFeCkKYGBgCmBgYHtyfQppbnN0YWxsLnBhY2thZ2VzKCJkZXZ0b29scyIpCmBgYAoKYGBge3J9CmxpYnJhcnkoZGV2dG9vbHMpCmluc3RhbGxfZ2l0aHViKCJXU3BpbGxlci9SYWRpYWxNUiIpCmBgYApgYGB7cn0KcGxvdChieS9ieCwgYngvYnlzZSwgeGxpbT1jKG1pbigoYnktMipieXNlKS9ieCksIG1heCgoYngrMipieXNlKS9ieCkpLCB5bGltPWMoMCwgbWF4KGJ4L2J5c2UpKSkKZm9yIChqIGluIDE6bGVuZ3RoKGJ4KSkgewpsaW5lcyhjKChieVtqXS0xLjk2KmJ5c2Vbal0pL2J4W2pdLCAoYnlbal0rMS45NipieXNlW2pdKS9ieFtqXSksCmMoYnhbal0vYnlzZVtqXSwgYnhbal0vYnlzZVtqXSkpCn0KYWJsaW5lKGg9MCwgbHdkPTEpOyBhYmxpbmUodj0wLCBsd2Q9MSkKYGBgCmBgYHtyfQpsaWJyYXJ5KFJhZGlhbE1SKQpyX2lucHV0IDwtIGZvcm1hdF9yYWRpYWwoYngsIGJ5LCBieHNlLCBieXNlKQpmdW5uZWxfcmFkaWFsKGVnZ2VyX3JhZGlhbChyX2lucHV0LCBhbHBoYSA9IC4wNSkpCmBgYAoKYGBge3J9CnBsb3RfcmFkaWFsKGl2d19yYWRpYWwocl9pbnB1dCwgYWxwaGEgPSAuMDUpKQpgYGAKCmBgYHtyfQpwbG90X3JhZGlhbChlZ2dlcl9yYWRpYWwocl9pbnB1dCwgYWxwaGEgPSAuMDUpKQpgYGAK